 classdef etsModel
    %etsModel Model for emissions trading
    %   Model for plant abatement with heterogenous pollution and emissions
    %   as a function of expenditure.
    
    properties
        
        % Data objects used in model estimation
        plants;
        bids;
        
        % Size of model
        S = 100;  % Number of simulations to form bidder expectations
        N = 156;  % Number of plants 
        T = 10;   % Number of compliance periods
        
        %% Model parameters
        
        %   Abatement function
        % costFunction = 'Cobb-Douglas';
        % alpha = [ 1 0.8 ]';
        % beta  = -1.1;
        
        costFunction = 'MaxEmissions';
        beta  = [ -1.6 -0.4 0.36 ];
        
        %   Firm observable distribution
        theta_h = [ 15 0.6 ];
        
        %   Time and firm effect distributions
        sigma_t  = 0.5; % Time effect distribution parameters
        sigma_i  = 1.0; % Plant effect distribution parameters
        sigma_it = 0.3; % Plant X period MAC shock std deviation
%         sigma_it   = 0; % Plant X period MAC shock std deviation
        
        %   Permit supply
        Qt;
        Pt;
        
        % Bidding parameters
        lambda    = 4;    % Parameter of Poisson distribution of bids
        sigma_k   = 0.30; % Standard deviation of log bid quantity 
                          %   relative to equilibrium bid quantity
        sigma_itk = 0.25 % Standard deviation in quantity prediction error
        
        % Command-and-control (emissions rate) parameters
        mu_ER;          % Mean of log emissions rate (initialized from cap)
        sigma_ER = 1.00 % Standard deviation of log emissions rate    
        
        %% Model endogenous objects
        period = 1; % To mark one period of interest, for subset / market clearing
        plant;
        bid;
        
    end
    
    properties (Dependent, SetAccess = private)
        plant_t; % Plant characteristics in a single compliance period
        bid_t;   % Plant X bid characteristics in a single compliance period
    end
    
    methods
        
        function obj = etsModel(  )
            % Initialize bidModel object
            
            % Store plant data
            obj.Qt = repmat(1000,1,10);
            obj.Pt = nan(size(obj.Qt));
            
            % Simulate plant characteristics 
            obj.plant = obj.drawPlants;
        end
        
        %% Retrieve plant and bid characteristics for a single period
        function plant_t = get.plant_t(obj)
            plant_t = obj.plant( obj.plant.id_period == obj.period, : );
        end
        
        function bid_t = get.bid_t(obj)
            bid_t = obj.bid( obj.bid.id_period == obj.period, : );
        end
               
        %% Sample plant characteristics
        function plant = drawPlants( obj )           
            str = RandStream('mt19937ar','Seed',20200820); % Random # stream

            % Draw plant characteristics (static)
            Hi  = exp(normrnd(obj.theta_h(1),obj.theta_h(2),obj.N,1));
            Ai_share = Hi / sum(Hi);  
            xii = normrnd(0,obj.sigma_t,obj.N,1);
            id_plant = [1:1:length(Hi)]';
            treatment = ones(length(Hi),1);
            plant = table(id_plant,Hi,Ai_share,xii,treatment);
            
            % Draw concentration limits C^bar_i for each plant
            plant.ERi = obj.drawPlantER( Hi );
            % KW: should we change to plant.ERit = obj.drawPlantER( Hi );?
            
            % Draw maximum emissions rate, which enters abatement cost
            plant.Ebari = obj.drawPlantEbar( Hi );
             
            % Expand from plant-level to plant X period level
            plant = plant( repelem(id_plant,obj.T), : );
            plant.id_period = repmat([1:obj.T]',obj.N,1);
            plant = movevars(plant,'id_period','After','id_plant');
            
            % Draw idiosyncratic plant X period pollution shocks eps_it
            plant.epsit = normrnd(0,obj.sigma_it,obj.N*obj.T,1);
            
            % Draw period shocks Xi_t
            xit = normrnd(0,obj.sigma_i,obj.T,1);
            plant.xit = xit( plant.id_period );
            
            % Assign permit allocation as a pro rata portion of 
            %   the period-specific emissions cap
            plant.Ait = plant.Ai_share .* obj.Qt( plant.id_period )';
            
            % Initialize emissions 
            plant.Eit = nan(size(plant,1),1);
        end
        
        % Draw sample permit bids
        %   Plants bid/offer a random number of bids at random quantities
        %   All {bid,qty} pairs are assumed to lie on plant-specific MAC curve
        function obj = drawBids( obj )
            obj = obj.initializeBids; 
            
            % Draw locations on marginal cost curve for each plant
            %   Log emissions centered on the equilibrium level
            str = RandStream('mt19937ar','Seed',20220304); % Random # stream
            
            bid_q_factor   = exp(normrnd(0,obj.sigma_k,size(obj.bid,1),1));
            obj.bid.bid_qE = obj.bid.Eit .* bid_q_factor;  % Emissions if bid executed
            obj.bid.bid_q  = obj.bid.bid_qE - obj.bid.Ait; % Bid quantity
            
            % Assign prices to each bid as equal to the marginal cost of
            %   abatement at that quantity plus an error term
            obj.bid.mac   = obj.solvePlantCostsMarginal( obj.bid.bid_qE );
            bid_p_factor  = exp(normrnd(0,obj.sigma_itk,size(obj.bid,1),1));
            obj.bid.bid_p = obj.bid.mac .* bid_p_factor;
            
            % Create dummy variable for bid in first half of compliance
            %   period, assumed equal to a constant
            obj.bid.D_first_half = ones(size(obj.bid.bid_p));
        end
        
        % Initialize plant X period X bid level data table
        function obj = initializeBids( obj )
            
            % Draw number of bids offered by each plant (Poisson)
            bid_count = poissrnd( obj.lambda, size(obj.plant,1), 1 );
            obj.plant = addvars(obj.plant,bid_count,...
                                'NewVariableNames','bid_count');
            
            % Initialize table at plant X period X bid level
            id_plant_period = [ 1:size(obj.plant,1) ]'; % Plant X period
            id_plant_period_long = repelem(id_plant_period,obj.plant.bid_count);
            
            % Create index for bid # within plant X period
            idx        = accumarray(id_plant_period_long(:),1,[],@(x){1:numel(x)});
            id_bid     = [idx{:}]';
            
            obj.bid = addvars(obj.plant(id_plant_period_long,:),id_bid,...
                'NewVariableNames','id_bid','After','id_period');    
            obj.bid = movevars(obj.bid,'bid_count','Before','id_bid');
        end
        
        % Draw plant concentration limits for command & control regime
        function [ ERi ] = drawPlantER( obj, Hi )
            
            % Emissions rate per plant
            ERbar       = (1/obj.N) * obj.Qt(obj.period) / mean(Hi); 
            obj.mu_ER = log(ERbar) - obj.sigma_ER^2/2;
            
            % Draw distribution of ER standards across plants
            ERi = lognrnd(obj.mu_ER,obj.sigma_ER,[obj.N,1]);   
        end
        
        % Draw maximum emissions for each plant
        function [ Ebari ] = drawPlantEbar( obj, Hi )
            
            % Maximum emissions at 10 * cap rate
            %   Note: this initialization implies that emissions cost 
            %     function will depend on the cap. In the data it will
            %     instead be fixed by an engineering formula and so
            %     invariant to the cap.
            Ebari = obj.Qt(1) * 10 * Hi / sum(Hi);
        end
        
        %% Emissions trading market equilibrium
        function [ obj ] = solveMarketAll( obj )
            % Solve for emissions market equilibrium in all periods
            for t = 1:obj.T
                
                % Update compliance period, which determines what period
                %   the market equilibrium will be solved *for*
                obj.period = t;
                
                % Market-clearing price
                obj.Pt(t) = obj.solveMarketOne;
                
                % Assign plant emissions
                obj = obj.assignPlantEmissions;
            end
        end
         
        function [ Pstar ] = solveMarketOne( obj )
            % Solve for emissions market equilibrium in one period given:
            %    - Plant characteristics (observable)
            %    - Plant characteristics (unobservable)
            %    - Market shock
            %    - Permit supply
            % return market-clearing price and aggregate emissions.
     
            % Parameters for market-clearing
            P0   = [ 0.001 1e6 ];
            Qcap = obj.Qt(obj.period);
            
            % Find market-clearing price (QD(P*) - QS = 0)
            objfun = @( P ) (obj.solvePlantEmissions( P ) - Qcap);
            Pstar  = fzero( objfun, P0 );     
        end
        
        function [ Et, Eit ] = solvePlantEmissions( obj, Pt )
            % Solve for plant-level emissions in a compliance period
            %   as a function of market prices.
            % Return: 
            %   Et  - Aggregate emissions 
            %   Eit - Vector of plant-period emissions
            
            switch obj.costFunction
                
                case 'Cobb-Douglas'
                                        
                    % Assign parameters
                    alpha0 = obj.alpha(1);
                    alphaH = obj.alpha(2);
                    betaZ  = obj.beta;

                    % Calculate emissions
                    C    = alpha0^(-betaZ/(1-betaZ))*(-betaZ)^(1/(1-betaZ));
                    E1i  = (obj.plant_t.Hi.^alphaH.*exp(obj.plant_t.xii)).^(1/(1-betaZ));
                    E2it = (exp(obj.plant_t.xit + obj.plant_t.epsit)).^(1/(1-betaZ));
                    EPit = Pt^(-1/(1-betaZ));

                    Eit = C * E1i .* E2it * EPit;
                    Et  = sum(Eit);
                    
                case 'MaxEmissions'
                    
                    % Assign parameters
                    beta0 = obj.beta(1);
                    beta1 = obj.beta(2);
                    beta2 = obj.beta(3);
                    
                    % Calculate emissions
                    E1   = exp( -( beta0 + obj.plant_t.xit + ...
                        obj.plant_t.xii + obj.plant_t.epsit ) );
                    E2   = obj.plant_t.Hi.^(-beta2);
                    EPit = Pt^(1/beta1);
                    
                    Eit = EPit .* [ E1 .* E2 ].^(1/beta1);
                    Et  = sum(Eit);
            end
            
        end
        
        function [ obj ] = assignPlantEmissions( obj )
            % Assign plant emissions in permit market equilibrium
            
            % Solve emissions for each plant
            [ ~, Eit_equil ] = obj.solvePlantEmissions( obj.Pt(obj.period) );
            
            % Assign emissions to the relevant compliance period
            tidx = (obj.plant.id_period == obj.period);
            obj.plant(tidx,:).Eit = Eit_equil;
        end
        
        %% Command & control status quo regulation
        function [ obj ] = solveCommandAll( obj )
            % Solve for plant-level emissions in all compliance periods
            %   as a function of aggregate emissions levels
            
        end
        
        function [ Eit ] = solveCommandPlantEmissions( obj, Et )
            % Solve for plant-level emissions in one compliance period
            %   as a function of the aggregate emissions level
            Et_baseline = sum( obj.plants.ERi .* obj.plants.Hi );
            Et_factor   = Et / Et_baseline;
            Eit         = Et_factor * (obj.plants.ERi .* obj.plants.Hi);
        end
        
        %% Cost and emissions functions used in both regimes
        function [ cit ] = solvePlantCostsMarginal( obj, Eitk )
            % Solve for plant-level costs as a function of bid quantities
            %   using each plant's marginal abatement cost curve
            % Eitk : Emissions at which to evaluate plant costs
            switch obj.costFunction
                
                case 'Cobb-Douglas'
                                               
                    % Assign parameters
                    alpha0 = obj.alpha(1);
                    alphaH = obj.alpha(2);
                    betaZ  = obj.beta;

                    % Calculate marginal costs at the implied emissions quantity
                    C    = alpha0^(-betaZ) * (-betaZ);
                    C1i  = ( obj.bid.Hi.^alphaH .* exp(obj.bid.xii) ).^(-betaZ);
                    C2it = ( exp(obj.bid.xit + obj.bid.epsit) ).^(-betaZ);
                    C3it = Eitk.^(betaZ - 1);

                    cit  = C .* C1i .* C2it .* C3it;
                    
                case 'MaxEmissions'
                    
                    % Assign parameters
                    beta0 = obj.beta(1);
                    beta1 = obj.beta(2);
                    beta2 = obj.beta(3);
                    
                    % Calculate marginal costs at the implied emissions qty
                    C1 = exp( beta0 + obj.bid.xit + obj.bid.xii + obj.bid.epsit );
                    C2 = Eitk.^beta1;
                    C3 = obj.bid.Hi.^beta2;
                    
                    cit = C1 .* C2 .* C3;
            end 
                
        end
        
        
        %% Functions to produce exhibits
        
        [ ] = plotAbatementAggregate( obj, filename, label, title );
        
        [ ] = plotEmissionsVsHeat( obj, filename, label, title );
        
        [ ] = plotAbatementBids( obj, plant_ids, filename, inlogs );

    end
    
 end
 
 
 
 
